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^ ABSTRACT 

o 

H 

\ We describe a search for close spectroscopic dwarf M star binaries using 

data from the Sloan Digital Sky Survey (SDSS) to address the question of the 
rate of occurrence of multiplicity in M dwarfs. We use a template fitting tech- 
nique to measure radial velocities from 145,888 individual spectra obtained for 
a magnitude-limited sample of 39,543 M dwarfs. Typically, the three or four 
spectra observed for each star are separated in time by less than four hours, but 
for ~ 17% of the stars, the individual observations span more than two days. 
In these cases we are sensitive to large amplitude radial velocity variations on 
time scales comparable to the separation between the observations. We use a 
control sample of objects having observations taken within a four hour period to 
make an empirical estimate of the underlying radial velocity error distribution 
and simulate our detection efficiency for a wide range of binary star systems. We 
find the frequency of binaries among the dwarf M stars with a < 0.4 AU to be 
3 — 4%. Comparison with other samples of binary stars demonstrates that the 
close binary fraction, like the total binary fraction, is an increasing function of 
primary mass. 



1. Introduction 



The study of the orbital motions of stars in multipl e systems has been an active area of 
research for more than a century (see iHearnshawl Il986l ). The fraction of stars which lie in 
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multiple systems (or conversely the fraction that are single), the dependence of these fractions 
on the mass of the primary, and the separation and mass ratio distributions of systems of 
multiple stars, are all observable quantities which can be measured through spectroscopic 
and imaging surveys and provide fundamental information in two important areas of stellar 
astrophysics. First, in certain circumstances observations of binary star systems directly 
measure the masses, radii, and luminosities of stars and provide the experimental bedrock 
for much of our theoretical understanding of the structure and evolution of stars. Second, 
surveys to determine the overall statistical properties of multiple systems within a population 
of stars help constrain their formation history; a detailed understanding of the fraction of 
binary stars across the HR diagram is an important part of a comprehensive picture of the 
process of star formation. To date, only a few measurements of the binary fraction and 
distribution for low-mass stars have been made - these stars are of low luminosity and the 
acquisition of the necessary large set of measurements is more difficult than for stars of higher 
mass and luminosity. They are of particular importance, however, because their formation 
and dynamical evolution within the systems in which they are born spans the mass range 
roughly between those of solar-type stars and planets. 

The current understanding of star formation is that stars form in small-N clusters which 
are t hen broken apart by (s hort-term) dynamical decay and (longer-term) dynamical destruc- 
tion (IGoodwin et al.ll2007l ). Due to the instability of multiple s ystems, some m embers of the 
system are likely to be ejected on a relatively short time scale; lAnosoval ( 119861 ) showed that 
in the vast majority of cases, the least massive star is ejected. On longer time scales , inter- 
actions with other stars in the star cluster disrupt binary stars (IGoodwin et al.ll2007| ). Both 
of these mechanisms predict that the multiplicity of stars should decrease to lower primary 
mass, but ejection is more important for close systems while interactions with other stars 
are more important for wide, loosely bound, systems. 

The properties of multiple star systems with Sun -like primaries have received much 
at tention of late. Recent work by iRaghavan et al.l (120101 ). which updates the seminal results 
of iDuquennoy &: Mayor! (Il99ll ). found that Sun-like stars have an overall multiple fraction 
of ~45%, with a log-normal distri bution in separati o n tha t peaks around 60 AU. A similar 
analysis of lower mass M stars by iFischer fc Marcyl (119921 ) found a somewhat lower overall 
binary frequency of 42±9%. For objects straddling the stellar-brown dwarf boundary, lAUen 
(200 7|) found that ul tracool dwarfs (M6 and later) have a binary frequency of 20±4%, while 
Sana fc Evansl ( 201QT) report that ~75% of O stars are binaries. Based on these results, Lada 
(120061 ) and IRaghavan et al.l (120101 ) point out that the overall fraction of single stars appears 
to be a decreasing function of stellar mass. 



Fischer fc Marcy 



(119921 ) report that 1.8±1.8% of early-M stars are binaries with 0.04 AU < 
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a < 0.4 AU based on a sample of 62 stars, while iBlake et al.l (1201 oh report that 2.5^' 6 



6/0 



of late-M and L dwarfs are binaries with a < 1 AU based on a sample of 43 objects. Ex- 



trapol ating from the overa l l sepa ration distribution functions presented by iRaghavan et al. 



(120 lOl ) and ISana &: Evansl (120101 ) provides estimates that 3.7% and 26% of G and O stars 
are binaries with a < 0.4 AU, respectively. A complete view of stellar multiplicity across a 
wide range of separations, however, requires the use of a variety of observational techniques. 
Nearby systems with wide separations may be directly resolved using high resolution imag- 
ing, while systems with small separations are most readily detected as spectroscopic binaries. 
In particular, measuring the clo se binary fractio n requ ires an extensive radial velocity (RV) 
survey of the type described by IRaghavan et al.l (120101 ) . 



Typically, individual multiple systems are identified in radial velocity surveys by fitting 
spectroscopic orbital solutions to velocity curves containing a large number of observations. 
Since close binaries are thought to be relatively rare for all but the most massive stars, a 
rigorous estimate of the multiplicity fraction requires a large number of observations of a 
large number of stars over long periods of time in order to detect an appreciable number of 
systems. This is particularly difficult for dwarf M stars because of their very low luminosities. 
An alternative approach is to use a small number of observations each for a very large number 
of stars to understand the rate of occurrence of multiple systems in a statistical sense (see 
Maxted fc Jeffries! 120051 ; iPourbaix et al.l 120051 ) . We take this approach in this paper, taking 
advantage of the enormous numbers of spectroscopic observations of M stars gathered during 
the course of the Sloan Digital Sky Survey (SDSS) to measure the close binary fraction 
of M stars using observations of 39,543 M dwarfs. While this data set is far larger than 
any used in previous research, it presents unique challenges. The data produced by the 
SDSS have significantly lower resolution and signal-to-noise ratio (S/N) then those typically 
used for making stellar velocity measurements. Additionally, SDSS typically only obtained 
three or four spectroscopic exposures of each object; in comparison. iFischer fc Marcyl (jl992j) 
obtained on average 15 observations per star. We demonstrate that an overall radial velocity 
precision of 4 km s -1 per observation can be achieved, more than sufficient for the detection 
of binary systems with short orbital periods, even with only few individual measurements. 
By quantifying the underlying statistical properties of the radial velocity measurements 
extracted from the SDSS spectra, and simulating the detection efficiency as a function of 
binary orbital separation and mass ratio, we make a robust measurement of the population 
of close binary M star systems. 
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2. SDSS Spectroscopic Data and Sample Selection 



2.1. The Sloan Digital Sky Survey 



The Sloan Digital Sky Survey (jYork et al.l |2000| ) made an imaging and spectroscopic 
survey of the sky using a large-f ormat CCD came ra (IGunn et al.l Il998l ) mounted on the 
Sloan Foundation 2.5 m telescope (IGunn et al.ll2006l ) at the Apache P oint Observatory New 
Mexico, to image the sky in five optical bands - u, g, r, i, and z ( iFukugita et al.lll996l ). 
The imaging data are reduced by a set of softwar e pipelines which produce a catalog of 



objects with cali brated magnitudes and positions (ILupton et al.l 12001 



Lupton et al.l 12003 



Hogget al.ll200l[ ) . Targets for spectr oscopy are selected from this catalog, mapped onto 
fiber plug plates (jBlanton et al.l 120031 ) . and observed with two dual fiber- fed spectrographs 
JUomoto et al.lll999h in a series of several 15-minute exposures. The spectra are optimally 
extracted, calibrated, combined (the combined spectra), and classified. 



2.2. Selection of the M Star Sample 



SDSS Data Release 7 (DR7: lAbazajian et al.ll2009l ) includes spectra of over 1.6 million 
objects, about 460,000 of which are stars. The parent M-star sample was selected from the 
spectroscopic data in DR7 by applying a series of cuts, first by magnitude and color (z < 19.5, 
i — z > 0.2 and r — i > 0.5). Next, data from plates with diffuse ionized gas emission and 
bad spectra (low S/N, missing data etc.) were removed, as were objects spectroscopically 
classified as stars earlier than K, as galaxies, or as quasars. The sample spectra were then 
examined by eye. Stars with obvious blue/white dwarf companions, stars superimposed on 
galaxies, and stars with the spectra of metal-poor (subdwarf M), K, or carbon stars were 
removed. The resulting sample (Knapp et al. in prepa r ation) contains 51,193 individual M0- 
L0 stars with spectral types as defined by lWest et al.l (120051 ). 



2.3. Multiple SDSS Spectra 



The SDSS usually obtained only one combined spectrum per object, but in a few cases 
there are two or more combined spectra, acquired one of two ways: (1) inadvertent observa- 
tions of the same object on more than one spectroscopic plate, or (2) duplicate observations 
for quality check purposes, either of the same object on two plates or on the same plate, re- 
plugged and re-observed. These observations can be used to search for varia bility over periods 
from days to years, an example being the search by iPourbaix et al.l (120051 ) for spectroscopic 
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binaries. In addition, the spectrum of each object is a combination of several exposures, and 
the spectra obtained from the individual exposures can also be used to search for variability, 
including radial velocity variability, as is done in the present paper. 

The SDSS spectroscopic data are acquired as a series of spectroscopic exposures observed 
sequentially with an exposure time of 15 minutes each (the individual 15 minute spectra). 
These spectra are then averaged to produce the combined spectrum. The data are taken this 
way both to increase the dynamic range of the spectroscopic observations and to allow for 
the easy removal of data artifacts such as "cosmic rays" . Because of the scientific value of the 
individual spectra, software to reduce and calibrate them was prepared for release in DR7 
along with the combined spectra. This data product is made available to enable studies of 
rapid spectral variability, since it provides, in almost all cases, three spectra observed with 
the same signal to noise ratio within a time period of one hour. The science enabled by 
these data includes identifying ob jects with short-peri od radial velocity variations such as 



compact degenerat e binary pairs flB adcncs et al 



lines in dMe stars flKruse et al.l I201CH : Hilton et al.ll2010h . 



2009[K and the variability of the emission 



The typical exposure sequence is three fifteen-minute observations within a time span 
of about an hour. However, in marginal weather more exposures may be required to achieve 
the required S/N, or the observational sequence may be interrupted by variable weather 
conditions or instrumental problems. In other cases, the required S/N may not be achieved 
in a single night, and the individual observations are separated by significantly longer times. 
This enables searches for spectroscopic variability in the time domain, ana l ogous to the 



photo metric variability studies th at have been carried out by lBlake et al.l (120081 ). iBhatti et al. 



( 120101 ) , and iBecker et al.l (120111 ) . The present paper uses the individual spectra to search 
for radial velocity variations in the large sample of M stars defined above. We use only 
the observations from the same plugging of a plate (see (3) above) and do not consider the 
repeat observations of individual stars or plates (cases 1 and 2 above) - these are discussed 
elsewhere. 



3. The Distribution of Radial Velocities 

3.1. The Final Samples 

We restricted the analysis to stars with 16 < i < 20.5. We broke the sample into sub- 
samples based on At, the total time baseline spanned by the observations of an object: a 
control sample with hr < At < 4 hr and an experimental sample with 2 d < At < 30 d. 
We expect very few stars in the control sample to undergo significant accelerations on such 
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short time scales. We use the statistical properties of the radial velocity measurements of 
the objects in this sample as an empirical estimate of the underlying radial velocity error 
distribution. Among the set of objects with 2 d < At < 30 d, the time spread between 
the observations is long enough that if the object is a close binary star we would observe 
significant radial velocity variations. 



3.2. Calculation of Radial Velocities 



We use a x 2 minimization technique to estimate rad ial velocities from the SPSS spec- 
tra employing the average low-mass star templates from iBochanski et al.l ( 120071 ) as a zero- 
velocity reference. Using a code written in Interactive Data Language (IDL), we determine 
the radial velocity that gives the best fit between the spectra and one of the templates by 
minimizing 

2 V^ r ^~ m ( A *)i2 



X 



D 



(1) 



where the sum is over all pixels in the spectrum (except for those flagged as bad, see below), 
fi is the calibrated flux of the i th pixel, m(Aj) is the value of the model at A, (the wavelength 
of the i th pixel), and a, is the estimated standard deviation of fi. We considered spectra 
only from the red arm of the SDSS spectrograph, which span s the region A = 5800 - 9200 A 
at a resolution of A/ A A ~ 1800 and contains 2 048 pixels (jStoughton et al.ll2002l ). Since 



the templates provided by IBochanski et al.l ( 120071 ) are normalized, there are two parameters 
that determine the model m that is the best fit for each spectrum: the radial velocity and 
an overall flux scaling. For each spectrum, we found the minimum of the x 2 curve using a 
iterative process in which we tested velocities in the range —1000 < RV < 1000 km s _1 at a 
resolution of 1 km s _1 and then on a finer grid about the velocity found to have the smallest 
X 2 in the previous step. This was repe ated twice with t he fin al step having a resolution 



of 0.1 m s _1 . The templates created by IBochanski et al.l (120071 ) give values for the flux at 
0.1 A intervals from 3825A to 9200A. For each test RV, the spectral temp late is interpolate d 
onto a Doppler-shifted wavelength grid using cubic spline interpolation (IPress et al.lll992l ). 
We fit each of the 11 templates (one for each spectral class from M0-L0) to each 15-minute 
spectrum of each star in our sample in order to find the best fitting template, and then select 
a single template that is used for all observations of a given target by finding the template 
that results in the lowest total x 2 when fit to all of the observations of a given object. We 
correct for the barycentric motion of the Earth by applying velocity corrections produced by 
the SDSS pipeline. 



During the fitting process, we exclude by down weighting spectral pixels that might 
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bias the resulting radial velocities. These include any pixels at wavelengths greater than 
A > 9150A, any pixels in the Ha emission region 6540 < A < 6585A, or pixels with the 
BADSKYCHI flag set by the SDSS pipeline. We exclude the reddest 50A of spectra since 
at these wavelengths telluric absorption and sky emission become significant and large radial 
velocity shifts may extend beyond the edges of the spectral templates. The Ha emission 
line can be very strong in active M stars and is known to vary s ignificantly even b etween 
individual observations of a given M star ( iBochanski et al.ll2007t iKruse et al.l 120101 ) . so we 
chose to exclude pixels within 22A of this feature. The BADSKYCHI flag indicates that 
sky emission lines are not being well fit by the SDSS spectral extraction pipeline, resulting 
in spectra that are potentially contaminated by sky emission. 



3.3. Identifying Binaries 

Given this sample of more than 145,000 radial velocity measurements, we hope to quan- 
tify our ability to detect short-timescale radial velocity variability in these data and estimate 
the rate of occurrence of short-period binary systems. Before doing this, we apply several 
cuts to the sample. We retain only radial velocity measurements from spectra that satisfy 
all of the following criteria: the average S/N of the pixels is greater than 10, the observation 
is not among the 10% of observations with the largest values of the \ 2 f° r the template 
fit, and after applying the two previous cuts, the spectra are of objects with at least three 
observations. We also remove objects from the sample that no longer fall into the definitions 
of our control and experimental samples due to the removal of observations by the previous 
cuts. 

After applying these cuts, 23,031 observations of 7,059 objects remain in the control 
sample and 6,845 observations of 1,452 objects in the experimental sample; thus, of the final 
set of objects analyzed, 17% have observations with a time span of two days or more. In 
Figured] we show the distributions of the number of observations of each object and At for 
both samples, as well as the distributions of z-band magnitudes and i — z colors. 



For each of the objects in both samples, we define relative velocities, ARVi 



AW - W £<HVi-(S/N)< (9) 

where RVj and (S/N)j represent the RV, after applying the barycentric correction, and the 
average signal-to-noise ratio of all pixels of the i th observation of the object, respectively. This 
is the difference between the radial velocity of the observation and the weighted average of 
the radial velocities of all observations where the weight function is (S/N)j. 
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Given only a small number of observations of an object, it is extremely important to 
understand the underlying radial velocity error distribution if we hope to reliably detect 
radial velocity variations. For example, systematic sources of error can result in significant 
non-Gaussian tails to the error distribution that could produce spurious detections of radial 
velocity variability. Instead of relying on statistical estimates of the radial velocity error for 
a given measurement, we use the distribution of ARV derived from the control sample as 
an empirical error distribution. 

We quantified the level of radial velocity variability for each object, x, as the sum of the 
absolute deviation of ARV 

where M is the number of observations of the object. We ran a Monte Carlo simulation of 10 7 
hypothetical sets of observations to determine a cutoff for x which is exceeded by only 10~ 3 
of the simulated objects. For each simulated object, we chose the number of observations 
from the distribution of the number of observations of each object in the control sample. 
We then randomly chose the ARV values for each simulated object from the 23,031 actual 
ARV values in the control sample and calculated x using Equation |3j From the results of 
this simulation, we determined a cutoff value of x > 10.4 km s _1 for variability significant at 
the 10 -3 level. Among the 1,452 objects in the experimental sample with 2 d < At < 30 d, 
22 exceed this cutoff and are therefore detected as radial velocity variables. These stars are 
listed in Table [TJ 

A histogram of the values of ARV computed for the control sample along with the 
best fit Gaussian distribution is shown in Figure [2j The standard deviation of the ARV 
values is 3.8 km s _1 , which compares favorably with t he rms velocity e r ror o f 5.5 km s _1 



at g — 18.5 and 12 km s _1 at g — 19.5 reported by lAbazajian et al.l ( 20091 ) . However, 
there is a small set of objects with very large ARV values. Among the 7,059 objects in the 
control sample, 14 have at least one observation for which ARV > 20 km s _1 . This is far 
greater than would be expected if the values of ARV followed a Gaussian distribution. We 
examined the spectra of these objects individually, but all appear to be normal M dwarfs. 
However, we did note that many of the spectra contained a small number of highly errant 
pixels, likely due to cosmic rays, that were not fully down weighted by the SDSS pipeline. 
We considered the possibility that these errant pixels were dramatically altering the radial 
velocity fit and to test this hypothesis, we refit the spectra of the 22 detected binaries in 
the experimental sample and the 14 outliers in the control sample with the 5 most deviant 
pixels removed from the fit. We found that the radial velocity changed by an average of only 
0.3 km s _1 , well within our radial velocity error, which allowed us to reject this explanation. 
Template mismatch also does not appear to be at fault as the average \ 2 °f these objects 
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differs only slightly from the mean for the whole sample. Based upon this, we conclude that 
the objects are either in fact undergoing very short-term radial velocity variability, possibly 
due to a very tight binary system, or that there is an error in the wavelength calibration 
for these spectra. The wavelength solutions for each spectrum taken with each SDSS fiber 
are determin ed through a combinati on of fits to sky emission lines and arc lamp lines. As 
described by lAbazajian et al.l ( 120091 ) . starting with DR7 these solutions were constrained 
to vary smoothly between fibers, resulting in overall calibration precision of ±2km s _1 and 
significantly reducing the rate of errant solutions. Still, it is possible that rare problems with 
the wavelength solutions may be responsible for a small number of our measured RV shifts. 
In any case, these objects are of interest and warrant further observation. They are listed 
in Table [2j We also note that a significant percentage of Sun-like stars are k nown to be in 
higher-order multiple systems (34% doubles, 9% triples; iRaghavan et al.ll2010l ). We are only 
sensitive to short period systems, so even in a triple system we would only detect the reflex 
motion of the tight inner pair. 

The raw numbers listed above show that 22 of the 1,452 stars (1.5%) in the experimental 
sample show statistically-significant radial velocity variations, compared to 14 of the 7,059 
(0.20%) of the stars in the control sample. If the small number of large radial velocity 
variations in the control sample is due to some unidentified instrumental or analysis problem, 
we would also expect 0.20% of the experimental sample (3 stars) to show these variations; 
thus the raw binary fraction in the experimental sample is measured at the la level (the 
on-average smaller number of observations available for the control sample is accounted for 
in the calculation of x in Equation [3]). 



3.4. The M Star Close Binary Fraction: a Bayesian Analysis 

The raw binary star fraction of 1.5% derived above underestimates the true binary 
fraction which could in principle be measured by radial velocities of the accuracy available 
from the SDSS spectra, due to eff ects of inclination an d the small number of observations 



of each star (see the discussion by lPourbaix et al.l 120051 ) . In this section, we use the radial 



velocity measurements of the stars in the control and experimental samples to estimate 
the true close binary fraction of the objects in our experimental sample by quantifying the 
number of objects that exhibit statistically significant radial velocity variations. To do this, 
we take a Bayesian approach to estimate the likelihood of a given close binary fraction given 
our observations and any prior information. 



Let D represent the number of radial velocity variables detected, Nc a value for the 
close binary fraction, and B knowledge of any relevant background information about the 
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objects and the observations of the objects, such as their magnitudes and the times of the 
observations. Bayes' theorem states that 

P(N C \D, B) cc P(D\N C , B) ■ P(N C \B) (4) 

where P(Nq\D, B), the probability of Nq given D and B, is the posterior distribution; 
P(D\Nc, B), the probability of D given Nq and B, is the likelihood distribution; and 
P(Nq\B), the probab ility of Nc given B, is the prior distribution for any Nc and D 



flSivia & Skillind 120061 ). 



As D represents the number of objects detected as radial velocity variables, it follows 

that 

P(D\N C , B) = J2 P({O h ,O l2 , . . . , O lD }\N c , B) (5) 

where P^O^, O i2 , . . . , O iD }\Nc, B) is the probability that exactly the D objects, i±, i 2 , ■ ■ ■ , %d-> 
are detected as radial velocity variables and the sum is over all sets of D objects. Since 
whether or not one object is detected is independent of whether or not another object is 
detected, 

P({O n , O i2 , . . . , O tD }\N c , B) = H P(0 3 \N c , B) ■ J] P(0 3 \N Cl B) (6) 

jeL jgL 

where L — {ii, i 2 , ■ ■ ■ , id}, P(Oj\N c , B) is the probability that the j th object is detected and 
P(Oj\N c , B) is the probability that the j th object is not detected. Substituting this into 
Equation [5] gives 

P(D\N C , 5) = 3[] P{0 3 \N C) B) ■ H PiOjlNc, B)} (7) 
jeL jgL 

We must now consi der how to calculate P(Oj\N c , B) and P(Oj\N c , B). Following 



Maxted fc Jeffries! (120051 ) in assuming that the only cause of radial velocity variability is 
stellar multiplicity, it follows that the probability that the j th object is detected as a radial 
velocity variable is NcPdetectj + (1 — Nc) ■ 10~ 3 where Pdetectj gives the probability that we 
will detect the j th object as a radial velocity variable if it is in fact a binary star. Note 
that the term (1 — Nc) ■ 10~ 3 arises from the fact that 1 — Nc is the probability that the 
object is not a binary star while 10~ 3 is the probability that an object will be detected 
as a radial velocity variable if it is not a binary star by the cutoff value of the variability 
metric we established to determine whether or not an object is a radial velocity variable. 
Thus, P(Oj\N c ,B) = N cPde tect,j + (1 - N c ) • 10~ 3 and P(0 3 \N C ,B) = 1 - P(Pj\N c ,B) = 
1 — [NcPdetectj + (1 — N c ) ■ 10~ 3 ]. We must now determine Pdetectj for each of the objects in 
the experimental sample, which can be done with a Monte Carlo simulation. 
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The radial velocity measurements described in section [3T31 and models for the distribu- 
tions of system orbital parameters formed the basis of a Monte Carlo simulation designed 
to determine Pdetectj, which is required to generate the posterior distribution. This Monte 
Carlo simulation consisted of 10 5 virtual binary stars for each of the 1,452 objects in the 
experimental sample, which enabled us to estimate the efficiency with which we detect bi- 
naries given a wide range of different binary systems. For this simulation, we drew binary 
parameters from distributions based on previous results in the literature. 

Semi-major axis, a: In order to have a reasonable chance to detect a binary in 
the experimental sample, th ere must be significa nt RV variations on the timescale of At, 
typically less than five days (IPourbaix et al.ll2005l ). Since the stars we are looking at are M 
stars, m 1 ~ m 2 < 0.5 M Q , and we are able to detect large RV variations, say ARV > 20 km 
s -1 , we are primarily sensitive to very close systems with a < 0.2 AU. Based upon our Monte 
Carlo simulations, we estimate that our average detection efficiency is 69% at a = 0.1 AU 
and 16% at a = 0.4 AU. The distribution of the semi-major axes for systems with such 
small separations is not well known. We consider two different distributions of a. The first 
is a uniform distribution from 0.01 to 0.4 AU, while the second is a linear distribution such 
that P{a) oc a that also runs from 0.01 to 0.4 AU. Note that we might expect the uniform 
distribution to overestimate the value of Pdetectj as the correct distribution should have 
a higher pro bability of larger separations and a corresponding lower probability of smaller 
separations ( lAllenll2007l ). This overestimate of p detect, j in turn implies that the best fit-binary 
fraction Nq will be underestimated. 



Mass ratio, q: We follow lAUenl ( 120071 ) in using a power law distribution with a mini- 
mum of q = 0.02. Thus, if we let 7 represent the power law index, the probability distribution 
function of q is 



P(q) 



(8) 



for .02 < q < 1 and P(q) 



I0.02 <? 7 d( l 

for < q < 0.02. We te st distributions with three different 



values for 7. We test 7 = 1.8, as found by|A 
values on the la confidence interval given by 



Alien! (|2007|). 



lenl (12007 ). and 7 = 1.2 and 2.2, the extreme 



Primary mass, m\. Unfortunately, there is not a good method for determining the 
mass of the primary from the SDSS spectra as there are no well-calibrated mass-color or mass- 
luminosity relationships using the SDSS filters. Additionally, not knowing the metallicity or 
age of the star increases the uncertainty in determining its mass. However, this is not a major 
issue because, while the optical spectra and brightness of cool M dwarfs cover a large range, 
the corresponding mass range is modest, and we are able to obtain a rudimentary estimate 
of the primary mass using relations from the literature. Using the color transformations 
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provided by lDavenport et al.l (120061 ). from the apparent magnitude of the star in the r and i 
ba nds we can estima te the i — J color. Also, using the color-magnitude relationships provided 
by I West et al.l (120051 ). we can determine the absolute magnitude, Mj from the i—z color of the 
ob jects. Since i — J — M j — Mj, we can estimate Mj. Finally, the mass-luminosity relations 



of iDelfosse et al.l (120001 ) allow for the estimation of the star's mass from Mj. While the 
individual values of the mass obtained by these means have large systematic uncertainties, 
this is mitigated by selecting primary mass from a uniform distribution from 0.75 to 1.25 
times the photometrically estimated mass in our Monte Carlo simulations. 

Eccentricity, e: It is known that binaries with very short periods (P<~10d) are highly 



likely to undergo tida l circularization (jDuquennoy &: Mayorlll99ll ; iMeibom fc Mathieull2005 



Raghavan et al.ll2010l ) and we can reasonably assume circular orbits, although the largest- 
separation pairs may have e ^ 0. 

Orbital phase, /: The orbital phase at the time of the first observation is chosen from 
a uniform distribution from to 2tt. 



Inclination, 

radians. 



The inclination is chosen from a uniform distribution from to 7r 



Longitude of periastron, u: The longitude of periastron is chosen from a uniform 
distribution from to 27r radians. 

With our Monte Carlo simulation we generate hypothetical radial velocity curves for 
binary systems given the actual times of the observations of each object and Keplerian orbits 
given randomly selected system parameters from the above distributions. We estimate errors 
on the individual data points in the simulated radial velocity curves by randomly selecting 
values from the ARV distribution of our control sample. We can then compute x using 
Equation |3l If the value of x is greater than the cutoff value determined in section |373| 10.4 
km s _1 , this system is tagged as detected. By running 10 5 trials for each of the 1,452 objects 
in our experimental sample and determining the fraction of trials in which we detect the 
simulated object, we estimate Pdetectj as a function of the binary system parameters. 

The prior distribution, P(Nq\ B), al l ows u s to incorporate any relevant outside informa- 
tion into our analysis. We follow lAUenl ( 120071 ) in choosing a prior that assumes no outside 
knowledge and therefore is not biased towards any particular v alues of Ng- Thus, as Nc is 



a scale parameter, the proper prior to use is the Jeffreys' prior (jSivia fc Skilling 120061 ) . 



P(N C \B) cc — . 



(9) 



which is equivalent to uniform in the log of N ( 



c- 



13 



Given the likelihood, P(D\Nc, B), and prior distribution, P(Nc\B), we can calculate 
the posterior distribution: 

P(N C \D, B) oc P(D\N C , B) ■ P(N C \B) (10) 
oc ^(11^1^-5) •IlWcB)] • -L (11) 



c 



YllTllNcPdeteetj + (1 - iV c )10- 3 ] ■ - {NcVdxtectj + (1 ~ A^HT 3 )]] ■ 



We can then determine the constant of proportionality using the normalization condition 

1 

P(N C \D, B) dN c = 1. (13) 

From this posterior distribution, we are able to determine the best fit value of the close 
binary fraction and a confidence interval on this value. 

Table [3] gives the average value of Pdetectj over all the objects and the best fit value 
of No and the la confidence interval on Nc for each combination of a and q distributions. 
From Table [3j it is clear that changing the power law index 7 has negligible impact on Nc- 
However, using a linear distribution for a as opposed to a uniform distribution increases 
Nq as expected, due to the decrease of Pdetectj that results from using a distribution with 
a higher probability of large separations. The increase in Nc that results from using a 
linear distribution instead of a uniform distribution is 0.9%. Figure [3] shows the posterior 
distributions for the close binary fraction for both a uniform and linear distribution of a. 
Assuming the uniform distribution in a, the close binary fraction of the objects in our 
sample is 2.9to'^%, while for a linear distributio n in a we find 3.8^o;g%. We note that since 



our sample is magnitude limited, Branch bias (IBranchlll976l ) may lead us to overestimate 
the rate of occurrence of multiple systems. For a population of binary systems with two 
identical stars, this can bias the measured binary fraction by up to 35%. We have not taken 
this into account in our analysis. 



Results and Discussion 



Previous work on the stati stical properties of mu ltiple star systems, such as t hat by 



puauennoy &: Mayor! (Il99ll ) and lRaghavan et al.l (120101 ). has largely excluded M stars. iFischer &: Marcy 
(119921 ) included spectroscopic binaries in their analysis of M star multiplicity and estimated 
that ~ 1.8% of M dwarfs are binaries with 0.04 AU < a < 0.4 AU. This result, given that it 
is based on a very small sample of targets, is fully consistent with our estimate of the binary 
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fraction across a similar range of separations. The rate of occurrence of multiple systems 
and the dependence of that rate on primary mass, may be an importa nt o bservational clue 



to a u nified theory of star formation. It has been shown by lLadal (120061 ) and lRaghavan et al. 



( 120 lOl ) that the overall multiple fraction is a strong function of stellar mass, decreasing from 
near 80% at 1OM to 20% for brown dwarfs with masses < O.1M . We gathered values from 
the literature for the binary and close-binary fractions as a function of primary mass. These 
are listed in Table HI 



For the overall multiple fr action, Nt, we ta k e the results of lAUenl (120071) for the brown 
dwarfs and lowes t mass stars, iFischer fc Marcvl (119921) fo r M st ars, iRaghavan et al.l (120101 ) 
for G stars, and Mason et al. ( 2009 ) and ISana k, Evans ( 2010 ) for O stars, for which we 
ass ume an average prima ry mass of 40M Q based on the studies of O star eclipsing binaries 
by IWeiden fc Vinki (120101 ). For the close binary fract i on, N g, of the lowest mass stars and 
brown dwarfs, we take the results from iBlake et al.l (120101 ) for the semi-major axis range 
a < 1 AU and scale down by a factor of 0.4 to estimate the binary fraction in the range 
a < 0.4 AU. For M stars we take the res ults presented he r e. Fo r G stars we integrate the 
overall semi-major axis distribution from IRaghavan et al.l (120101 ). inferred from the period 
distribution by assuming binary systems with a total mass of 1.5M , for a < 0.4AU. For B 
stars we take the lower limit on the rate of occurrence of spectroscopic binaries with periods 
less the lOd derived bv lMazeh et al. (12 006) by applying the most extreme correction for the 
Branch bias to the result of Wolff ll978h and assume a stellar mass of 3.8M . Finally, for 
O stars we use the empirical Cumulative Distribution Function (CDF) of observed O star 
spectroscopic binary periods to estimate Ng for a < 0.4AU assuming the binaries have a 
total mass of 70M Q , and therefore periods less than lOd. The relation between Nt and Nc 
and primary mass is shown in Figure HI 

These results show that the binary fraction is lower for lower mass stars across a wide 
range of orbital separations. However, the ratio Nc/N? appears to increase with increasing 
primary mass, indicating that either the peak of the semi-major axis distribution moves 
to smaller values of a for larger stellar masses or that the overall shape of the separation 
distribution evolves significantly with primary mass. Based on well-studied samples of M/L 
and G stars , there is evidence that in fact the lowest mass binaries ten d to have smaller 
separations (jDuquennoy &: Mayorlll99ll ; lAllenll2007l ; IRaghavan et al.ll2010l ). We caution that 
for O and B stars the multiplicity studies are not as observationally complete as those 
focusing on L and G stars, though with such high multiple fractions for these massive stars 
already reported in the literature it is unlikely that a sig nificant p o pulat ion of wide separation 
systems remains undetected. Even if, as suggested by iKroupal (120111 ). there is a universal 
Initial Period Function (IPF) for all binary stars, the processes of dynamical destruction by 
stars external to bound stellar systems and dynamical decay, i.e. ejection by stars internal to 



the s t ellar system , must play an important role in the evolution of binary properties ( jHeggie 
19751 . lHillslll975l ) and a wide range of binary properties for field stars may be expected. 



5. Conclusions 

Using a large set of spectroscopic observations of cool dwarf M stars from the SDSS 
that provide sparsely-sampled time series (typically 3-4 observations) we isolate 22 close 
spectroscopic binary candidates by radial velocity variability. Comparison with the radial 
velocity distributions of a control sample also from SDSS shows that we have detected the 
raw binary fraction at a 7a level. The cadence of the observations, and the SDSS radial 
velocity accuracy, make us sensitive to close binaries, with separations less than about 0.4 
AU. The detectability of spectroscopic binaries in the SDSS observations is evaluated using 
a large Monte Carlo simulation of the observational properties of the underlying population. 
Taking into account this detection efficiency and the total number of objects in our sample, 
the close binary fraction (a < 0.4 AU) of dwarf M stars is determined to be 2.9+q g% assuming 
a uniform prior on the distribution of orbital separation a. Our estimated close binary 
fraction of 3 — 4%, depending on the prior for a, is broadly consistent with earlier results on 
the statistical properties of low-mass multiple systems. 

By comparing our measurement of the close binary fraction of M dwarfs to previous 
results in the literature, we have shown that the close binary fraction, like the overall binary 
fraction, is an increasing function of primary mass. This result has implications for the 
formation of, and the early dynamical evolution of, systems of low mass stars and indicates 
that the overall distribution of binary separations may be a strong function of stellar mass. 
In the future, the methods described in this work can be used to investigate the binary 
fraction in the whole sample of stars observed spectroscopically by SDSS - over 460,000 
stars, mostly of types K, G and F, and members of the thick disk and halo. 
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Fig. 1. — : Subfigures (a)-(f): histograms of the number of observations of each object in the sample with 
hr < At < 4 hr, the number of observations of each object in the sample with 2 d < At < 30 d, 
the i magnitude in the combined sample, the i — z color in the combined sample, At in the sample with 
hr < At < 4 hr, and At in the sample with 2 d < At < 30 d, respectively. 
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Table 1. Object name, i magnitude, value of x calculated using Equation El the range in 
RV values (RV range =RV max — RV min ), and At for each object in the sample with 
2 d < At < 30 d that we detected as a binary. 



Object Name 


i 


X 


FOVr (in ge 


At 






(km/s) 


(km/s) 


(d) 


SDSSJ002228. 76-010806.9 


16.59 


12 


38 


19.002 


SDSSJ004747.84-004550.5 


17.02 


17 


38 


17.938 


SDSSJ010122.67-005311.6 


17.07 


11 


35 


2.019 


SDSSJ011038.87-011409.5 


17.31 


12 


32 


6.013 


SDSSJ011054.12+005227.4 


16.54 


14 


37 


6.013 


SDSSJ074428.29+191554.7 


17.47 


70 


251 


5.035 


SDSS J074646. 12+282629.9 


17.55 


12 


41 


7.056 


SDSSJ083027.95+454736.5 


17.33 


13 


35 


4.184 


SDSSJ083723.40+141211.1 


16.95 


11 


26 


14.936 


SDSSJ084841. 17+232051.7 


16.34 


92 


202 


8.090 


SDSSJ105030.21+421451.4 


16.44 


24 


90 


3.004 


SDSSJ105715. 76+430945.9 


16.33 


24 


97 


3.015 


SDSSJ114030.06+154231.5 


16.18 


47 


98 


17.963 


SDSSJ114050.85+532304.0 


16.65 


13 


35 


24.966 


SDSSJ115124.34+371953.8 


16.37 


14 


40 


2.135 


SDSSJ121944.13+260759.8 


16.13 


12 


36 


2.865 


SDSSJ163215. 69+005918.6 


16.66 


24 


98 


24.994 


SDSSJ163401. 41+005010.0 


16.87 


14 


53 


24.994 


SDSSJ204845.85+004001.3 


17.47 


47 


182 


19.957 


SDSSJ212546.00-060858.6 


16.12 


16 


43 


22.955 


SDSSJ220848. 44+000409.9 


17.12 


23 


59 


20.869 


SDSSJ225450.30-101003.2 


17.17 


17 


10 


21.975 



Table 2. Object name, i magnitude, value of x calculated using Equation EJ the range in 
RV values (RV range =RV max — RV min ), and At for each object in the sample with 
hr < At < 4 hr with at least one observation for which ARV > 20 km/s. 



Object Name 


i 


X 


RVrange 


At 








(km/s) 


(km/s) 


(d) 


SDSSJ003144.47- 


h003033.6 


17.49 


14 


40 


0.758 


SDSSJ021838.39- 


h004946.8 


16.98 


16 


41 


0.685 


SDSSJ072642.62- 


K414243.2 


16.64 


14 


44 


1.188 


SDSSJ075243.70- 


K254928.3 


17.39 


15 


41 


0.723 


SDSSJ082833.58- 


h341531.6 


16.79 


15 


35 


0.554 


SDSSJ085921.75- 


h371147.9 


18.19 


15 


34 


1.087 


SDSSJ092345.54- 


^222432. 4 


16.56 


18 


133 


1.242 


SDSSJ111647.81- 


^294602. 7 


16.24 


14 


36 


0.585 


SDSSJ125508.85- 


H320849.9 


17.04 


10 


33 


1.053 


SDSSJ143524.64- 


K232249.5 


16.97 


17 


48 


0.564 


SDSSJ154848.35- 


h 362803. 7 


17.04 


17 


42 


0.438 


SDSSJ163020.19- 


h305254.5 


16.62 


20 


55 


1.343 


SDSSJ163355.96- 


h293725.0 


16.98 


13 


34 


1.343 


SDSSJ163544.45- 


h243032.3 


16.78 


13 


35 


1.873 



Table 3. ip detect,]) i the best fit value of Nc and la (68.3%) confidence interval on Nc for 

each combination of a and q distributions. 



P(a) P(q) (p d etect,j ) N c 



Uniform 


7 = 


1.2 


0.45 


3.ol?:jj ( . 


Vc 


Uniform 


7 = 


1.8 


0.46 


o q+o.6< 


% 


Uniform 


7 = 


2.2 


0.47 


2 q+°- 5 < 


% 


Linear 


7 = 


1.2 


0.34 


4.0i°;« ( 


% 


Linear 


7 = 


1.8 


0.35 


q o+0.9( 
°-°-0.9- 


Vo 


Linear 


7 = 


2.2 


0.36 


q o+0.8( 


Vo 



Table 4. Summary of current results for the total binary fraction and close binary fraction 
for stars of various primary masses. See Section H] for information on the sources of these 

values. 



Primary Mass 


N T 


N c 


(M ) 


(%) 


(%) 


0.1 


20 


1.0 


0.4 


42 


2.9 


1.0 


45 


3.7 


3.8 




>8.0 


40.0 


75 


26.0 
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Fig. 2. — : Histograms of ARV and ARV/a for the control sample along with the best fit Gaussian distri- 
butions. 




A RV (km/s) ARV/a 



Fig. 3. — : Posterior distributions for a uniform (solid line) and linear (dashed line) distribution of a, 
respectively. 
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Fig. 4. — : Values of the close binary fraction and total binary fraction given in Section 0] and shown in 
tabular format in Table [4] along with best fit line for each fraction and total binary fraction. 
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